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AN  APPROXIMATE  FAST  FOURIER  TRANSFORM  TECHNIQUE 
FOR  VERNIER  SPECTRAL  ANALYSIS 

INTRODUCTION 


To  detect  the  presence  of  very  narrowband  weak  signals  in  noise,  and  to 
measure  their  center  frequencies  accurately,  it  is  necessary  to  Fourier  trans¬ 
form  a  long  time  segment  of  the  available  process.  When  the  center  frequencies 
of  the  signal  components  are  unknown  and  the  total  search  bandwidth  of  interest 
is  large,  this  procedure  demands  storage  and  computation  of  many  degrees  of 
freedom,  that  is,  search  of  a  large  time -bandwidth  product  space.  It  would  be 
advantageous  if  a  quick,  coarse  search  for  narrowband  components  could  be 
conducted,  followed  by  a  finer  vernier  analysis  over  a  limited  band  where  the 
presence  of  narrowband  components  has  been  indicated.  Such  an  adaptive  pro¬ 
cedure  would  be  less  time-consuming  and  require  less  storage.  Also,  if  the 
procedure  did  not  need  to  be  exact,  but  yielded  an  approximation  with  accept¬ 
able  sidelobes,  the  required  storage  and  computation  might  be  reduced  further. 

This  report  presents  just  such  a  technique,  which 

1.  accepts  the  input  process  in  smaller  time  segments  as  they  are  available, 

2.  performs  a  reasonable-size  weighted  fast  Fourier  transform  (FFT)  on 
each  overlapped  segment, 

3.  stores  only  that  frequency  portion  (at  each  segment)  where  narrowband 
components  are  indicated  to  be  present,  and 

4.  performs  a  small-size  weighted  FFT  over  the  total  data  record  available, 
for  each  frequency  bin  stored. 

Steps  1  and  2  permit  smaller-size  FFTs  than  would  be  required  if  the  total  data 
record  were  spectrally  analyzed  in  one  operation.  Steps  3  and  4  constitute  the 
adaptive  feature  of  this  technique.  The  last  transform  over  time  (delay)  in  step  4, 
for  each  frequency  bin,  is  a  vernier  frequency  analysis,  measured  from  the 
center  of  each  bin;  the  degree  of  approximation  of  this  technique  is  the  subject 
of  this  report. 
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Some  past  work onperforming  large-size  FFTs  by  means  of  several  smaller 
FFTs  is  reported  in  references  1  and  2.  The  methods  reported  there  are  exact, 
but  they  consume  more  time  and  require  more  storage  than  the  method  to  be 
presented  here.  In  particular,  the  two  methods  of  reference  1  require  too  many 
small-size  FFTs,  and  the  method  of  reference  2  requires  additional  multiplica¬ 
tions  by  complex  exponentials  and  a  fair  amount  of  storage.  The  approximate 
technique  of  reference  3  is  similar  to  the  one  outlined  above,  up  to  step  4,  with 
the  notable  exception  of  overlapped  weighting;  at  that  point  the  technique  of  ref¬ 
erence  3  requires  transformation  back  to  the  time  domain  followed  by  another 
transform  to  the  desired  frequency  domain.  Additional  transforms  are  required 
in  this  last  technique,  and  it  produces  greater  sidelobes  than  the  new  technique, 
especially  when  the  temporal  weighting  is  judiciously  selected. 

FUNDAMENTAL  SPECTRA L  RELATIONSHIPS 
LARGE-SIZE  FFT  APPROACH 

Before  embarkiiig  on  the  approximate  technique,  we  review  the  standard 
large-size  FFT  approach  to  spectral  analysis.  Suppose  a  data  waveform  x(t) 
is  sampled  at  time  instants  nA,  n  integer.  Then  the  voltage  density  spectrum 
that  can  be  computed  is* 

V(f)  s  f  dt  exp  (-i2irft)  x(t)  u(t)  A  64(t) 

4  (i) 

=  A  £  exp  (-i2irfnA)  x(nA)  u(nA)  , 
n 

where  u(t)  is  a  temporal  weighting  deliberately  imposed  to  control  spectral 
sidelobes,  as  will  be  discussed  shortly;  see  figure  1A.  The  finite  duration  of 
u(t)  terminates  the  integral  and  sum  in  equation  (1)  at  finite  limits.  The  im¬ 
pulse-train  function  in  (1)  is  defined  as  +he  infinite  sum 

«A  (t)  =  £  6  (t  -  nA)  .  (2) 

n 

The  integral  representation  in  (1)  allows  us  to  expresst 

V(i)  =  X(f)  •  U(f)  •  $£<9  ,  (3) 


♦All  integrals  are  over  the  range  of  nonzero  integrand. 

+The  Fourier  transform  of  the  lower-case  time  function  x(t)  is  the  upper¬ 
case  frequency  function  X(f);  this  notation  is  used  throughout. 
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where  •  denotes  convolution,  and,  in  keeping  with  (2), 

£,(f)=  E  S(f~).  «> 

m 

Thus,  the  observed  spectrum  V(f)  is  the  convolution  of  X(f)  with  the  set  of 
windows*  U(f)  •  6y4(f),  which  is  depicted  in  figure  IB. 


Because  the  ideal  spectral  window  is  a  single  impulse  at  f  =  0,  the  aliased 
mainlobes  at  m/A,  m  ^  0,  are  undesired.  Also,  the  window  U(f)  is  desired 
to  be  narrow,  with  very  small  sidelobes.  Since  the  weighting  u(t)  is  of  limited 
duration,  the  mainlobe  width  of  U(f)  is  not  zero,  but  is  inversely  proportional 
to  the  time  duration. 

Now,  if  the  voltage  density  spectrum  V(f)  is  computed  at  multiples  of 
(LA)”*,  where  L  A  is  the  time  duration  of  u(t),  we  obtain 

V  (— 1=  A  Y\  exp  (-i2irnp/L)  x(nA)  u(nA)  ,  p  integer  .  (5) 

\LA/  q 


*In  the  time  domain,  u(t)  is  called  a  weighting;  the  corresponding  Fourier 
transform  in  the  frequency  domain,  U(f),  is  called  a  window.  This  nomencla¬ 
ture  is  used  throughout. 


A 
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Since  V(f)  is  periodic  of  period  1/A  (see  (1)),  (5)  need  be  computed  at  L 
different  points;  thus  it  can  be  realized  as  an  L-point  FFT  of  sequence 
jx(nA)  u(nA)|.  For  line  frequency  analysis  (that  is,  large  LA)  the  size  L  of 
the  FFT  may  be  too  large  to  compute  easily,  under  storage  and  time  limitations. 
The  values  in  (5)  are  samples  of  tae  convolution  of  figure  IB  with  voltage  density 
spectrum  X(f). 

APPROXIMATE  FFT  TECHNIQUE 


Just  as  we  started  above  with  an  integral  definition  of  a  spectrum,  then 
showed  that  samples  of  this  spectrum  were  attainable  with  an  FFT,  we  begin 
with  the  spectral-delay  function,  a,  defined  as 

a(f,r)  =  f  dtexp  (— 12 «•  ft)  x(t)  w(t  -  r)  (t) 

(6) 

=  A  53  exP  (~12»fnA)  x(nA)  w(nA  -  r)  . 
n 


The  temporal  weighting  w  is  now  delayed  by  r  seconds;  if  the  duration  of  w 
is  Lw  seconds,  the  function  w(t  -  r)  picks  out  a  delayed  portion  of  data  x 
of  length  Lyy,  and  subjects  it  to  the  same  transform  as  in  (1).  This  operation 
is  depicted  in  figure  2A,  where  the  temporal  weighting  can  be  located  at  aa, 
bb, . . . ,  cc.  This  figure  is  drawn  for  50  percent  overlap  of  the  temporal  weight¬ 
ings;  however,  other  overlaps  are  possible  and  recommended  in  some  cases. 


The  next  step,  consistent  with  step  4  in  the  Introduction,  is  to  perform  a 
Fourier  transformation  on  the  delay  variable  r,  while  holding  frequency  vari¬ 
able  f  fixed.  The  general  definition  is  the  vernier  spectrum 


Y(f,  «>)  =  j  dr  exp  (-12  »  v  r)  a(f,  r)  d(r)  S65(t) 

=  S  53  exP  (-i2*pkS)  a(f,  kS)  d(kS)  , 
k 


(7) 


where  v  is  the  vernier  frequency,  d(r)  is  called  the  delay  weighting,  and 
S  is  the  separation  increment  in  delay  r  at  which  a(f,  r)  must  be  computed; 
that  is,  S  is  the  shift  between  temporal  weighting  locations  (see  figure  2A). 
Since  the  separation  S  in  delays  can  be  taken  to  be  smaller  than  the  temporal 
weighting  duration  ,  (7)  allows  for  overlapped  weighted  transformation  of 
the  available  data  (see  (6)). 
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The  operation  described  by  (7)  is  depicted  in  figure  2B.  When  the  Fourier 
transform  (6)  on  weighted  time  segment  aa  is  completed,  the  set  of  frequency 
components  denoted  by  the  vertical  line  of  Xs  at  the  end  of  segment  aa  are 
available.  Similarly,  frequency  component  values  at  the  ends  of  segments  bb, . . . , 
cc  are  indicated  in  figure  2B;  these  components  correspond  to  delayed  locations 
of  the  temporal  weighting  w.  Now,  for  a  fixed  frequency,  say  f^ ,  the  array 
of  (delayed)  frequency  components  indicated  in  a  horizontal  box  in  figure  2B  is 
subjected  to  a  delay  weighting  and  is  Fourier  transformed  according  to  (7), 
thereby  yielding  vernier  spectrum  Y(f},  »).  Similar  outputs  are  available  for 
other  (actyacent)  frequencies  of  interest,  such  as  f2  or  f3  . 


<-w  *4 


Figure  2A.  Overlapped  Temporal  Weightings 


FriqiMncy 


Figure  2B.  Adjacent  Delay  Weightings 


Figure  2.  Temporal  and  Delay  Weightings 
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Samples  of  the  vernier  spectrum  in 
is  the  duration  of  delay  weighting  d(r), 


9  at  multiples  of  (MS)"1, 
are  given  by 


where  MS 


Y  exp  (-i2*mk/M)  a(i,  kS)  d(kS), 
k 


m  integer, 


(8) 


whl  ih  can  be  realized  as  an  M-point  FFT  of  the  sequence  ja(f,  kS)  d(kS)|  of 
length  M.  The  periodicity  of  Y(f,  p)  in  p  ,  of  period  l/s  (see  (7)),  means 
that  (8)  need  be  computed  only  at  M  different  values  of  m. 


Finally,  samples  of  vernier  spectrum  Y  in  frequency  f  at  multiples  of 
(NA)-1  are  given  by  (using  (8)) 


/  P  m  \  /  P  \ 

r( - , - 1  =  S  Y  exP  (-i2»mk/M)  al —  ,  kSl 

\Na  MS/  k=0  \NA  / 


d(kS)  , 


(9) 


p  0, 1, « •  • ,  N  -  li  m  0, 1, « « « ,  M  ~  1  , 


where  delay  weighting  d(r)  has  been  selected  so  that  samples  |d(kS)|  are 
nonzero  only  for  k  =  0, 1, . . . ,  M  -  1.  The  values  of  a  needed  in  (9)  are  (using 
(6))  given  by 


Y  exp  (-i2«rpn/N)  x  (nA)  w(nA 
n 


-  kS)  . 


(10) 


In  order  to  put  (10)  directly  in  the  form  of  a  standard  FFT,  we  assume  that  the 
delay  separation  S  is  taken  as  an  integer  multiple  of  the  sampling  increment  A: 


S  =  IB  A  .  (11) 

Then,  if  temporal  weighting  w  has  nonzero  samples  |w(nA)f  onlyfor0<n<  N-l, 
(10)  becomes 


N-l 

exp(-i2»pkIg/N)  A  £  exp  (-i2*pm/N)  x(mA  +  klg  A)  w(mA), 
m=0 

(12) 


0  <  p  <  N  -  1,  0<k<M-l  . 
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The  exponential  phase  factor  preceding  the  sum  (FFT)  in  (12)  takes  on  particu¬ 
larly  simple  forms  for  two  special  cases  of  delay  separation  S:  For  Is  =  N/?, 
delay  S  is  equal  to  half  the  temporal  weighting  duration  Lw  and  is  termed 
50  percent  overlap;  for  I  =  N/4,  delay  S  is  one-quarter  of  1^  and  is  termed 
75  percent  overlap.  For  these  two  cases, 

50  percent  overlap,  I8  =  N/2,  phase  factor  =  (-l)pk  ;  (13A) 

75  percent  overlap,  Is  =  N/4,  phase  factor  =  (-i)pk  .  (13B) 

By  proper  branching  in  a  computer  program,  no  storage  or  complex  multiplica¬ 
tions  are  necessary  to  incorporate  these  phase  factors  in  (12),  prior  to  its 
usage  in  (9).  (An  alternative  approach  that  completely  circumvents  the  phase 
factor  in  (12)  is  described  in  appendix  A.) 

Equations  (12)  and  (9)  are  the  essential  results  of  interest.  We  now  inter¬ 
pret  thembymean8of  simple  examples  that  will  enable  us  to  make  good  choices 
of  temporal  weighting  w,  delay  weighting  d,  and  separation  (overlap)  S. 

INTERPRETATION  OF  THE  VERNIER  SPECTRUM 

In  appendix  B,  the  vernier  spectrum  is  shown  to  be  given  in  terms  of  X  by 


where  all  the  convolutions  are  on  r,  with  f  held  fixed.  D(v)  is  the  delay 
window  corresponding  to  the  delay  weighting  d(r). 

The  linearity  if  the  two  Fourier  transforms,  (6)  and  (7),  on  the  data  x(t) 
indicates  that  we  can  investigate  the  behavior  for  data  components  separately 
and  merely  add  the  results.  The  fundamental  component  is 

x(t)  =  exp  (i2wf0t),  X(f)  =  *(f  -  fo)  .  (15) 

At  this  point,  we  shall  make  a  series  of  reasonable  assumptions  and  require¬ 
ment.  ,  and  deduce  desirable  prope  '•ties  about  the  weightings  and  separations. 
The  first  assumptions  are 

(a)  excitation  frequency  fQ  <  (2  A)"1, 

(b)  coarse  analysis  frequency  f  <  (2A)"1, 

(c)  Lyy  ^  A  • 
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Assumption  (a)  av>  da  aliasing,  (b)  restricts  analysis  to  the  fundamental  range, 
and  (c)  requires  tin.  temporal  w  eighting  to  cover  many  samples  of  the  process 
x(t) .  Furthermore,  if 

(d)  temporal  window  W  has  low  si  delobes, 

tiie  only  term  in  the  sum  in  (14),  after  substituting  (15),  that  has  substantial 
value  is  that  for  m  -  0,  and  it  yields 

Y(f»  p)  3  [w(f  -  fo)  D(a  +  f  -  f  .  (16) 

A  plot  of  this  equation  versus  vernier  frequency  y  Is  given  in  figure  3,  where 
is  the  length  (duration)  of  delay  weighting  d(r ).  The  narrow  lobe  at  p  =  fQ  -  f 
is  the  desired  component;  this  component  in  displaced  from  the  coarse  analysis 
frequency  f  (corresponding  to  p  =  0)  by  -  f  Hz,  which  places  it  at  absolute 
frequency  f  +  (^  -  f)  =  fQ,  as  desired.  The  shape  of  this  lobe  is  governed  by 
the  delay  window  D;  thus,  if 

(e)  Ld»  8, 

(f)  delay  window  D  has  low  sidelobes, 

the  large  lobes  separated  by  1/S  Hz  in  figure  3  will  not  overlap  significantly, 
and  potentially  good  spectral  estimation  is  possible.  The  necessity  of  delay 
weighting  is  made  obvious  by  these  observations. 

A  Ylf.J/l 


Figure  3.  Vernier  Spectrum 
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There  are  a  few  additional  points  worth  noting  about  figure  3.  The  peak  height 
of  the  lobes,  W(f  -  fj  D(0),  is  a  function  of  the  exact  location  of  the  excitation 
frequency  fo  and  the  coarse  analysis  frequency  f.  This  undesirable  picket- 
fence  effect’*1  (which  was  not  present  in  figure  1)  can  be  minimized: 

(g)  choose  analysis  frequencies  jf^}  closely  spaced  (see  figure  2B). 


Then  Iffc-fol  is  small  for  some  value  of  k.  Also,  since  the  width  of  the  lobes 
in  figure  3  is  l/l^j,  where  is  the  delay  weighting  duration  and  will  betaken 
of  the  order  of  the  total  record  length  available  or  utilized,  very  fine  resolution 
in  p  is  possible.  Hence,  narrowband  components  closer  than  l/Lyy,  the  resolu¬ 
tion  capability  of  a  single  time  segment,  can  be  resolved  by  using  this  technique. 

FFT  CONSIDERATIONS 


Samples  of  the  vernier  spectrum  Y(f,  i>)  were  given  in  (9).  The  locations 
of  the  samples  aret 


f: 


r: 


i  2  N-l 

n  . . 

Vf  —  f  —  f  •  •  •  9  -  9 

NA  NA  NA 

1  2  M-l 

0,  ,  —  ,  •  •  • ,  — —  • 

MS  MS  MS 


(17) 


The  range  covered  by  the  vernier  frequency  v  is  S-1,  and  will  be  greater  than 
the  increment  in  fs  which  is  (NA)-*,  if  overlapped  temporal  weighting  is  used. 
And  since  the  full  range,  S-*,  would  encompass  a  spurious  lobe  for  values  of 
|^  -  f|  near  (2S)-1  (see  figure  3),  overlapping  is  necessary. 

The  approach  adopted  here  is  to  utilize  all  the  samples  in  f  at  separations 
of  (NA)-1,  and  use  only  samples  in  v  which  cover  a  range  of  (NA)-*;  that  is, 
we  use  the  central  portion  of  Y  cente  ed  around  r  =  0,  including  negative  fre¬ 
quencies.  In  terms  of  figure  2B,  adja<  ent  delay  weightings  at  fj,  f2,  f3  will  be 
employed.  The  alternative  time-savinj  procedure  of  attempting  to  utilize  all  of 
the  M  samples  in  v,  and  using  only  enough  samples  in  f  to  fill  in  the  frequency 
axis,  can  lead  to  a  very  bad  picket-fence  effect,  in  addition  to  large  spurious 
lobes  at  undesired  frequency  locations.  These  conclusions  follow  upon  piecing 
together  several  vernier  spectra  like  figure  3  for  appropriate  values  of  f  and 
excitation  frequency  fQ . 


♦See  reference  4,  page  47. 

tThe  upper  half  of  the  array  of  nu  nbers  in  (17)  corresponds  to  negative  fre¬ 
quencies.  Thus  thelast  samples  in  each  array  correspond  to  f  =-(NA)_1  and 
v  -  -(MS)”*,  respectively. 


\ 
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EXAMPLES 

The  general  guidelines  furnished  in  the  previous  section  do  not  yec  enable 
us  to  make  quantitative  selection  of  good  weightings  for  different  degrees  of 
overlap.  To  make  this  selection,  several  examples  are  considered  and  compared. 
The  numerical  examples  utilize 

1  1 

A  =  - seconds,  N  =  1024,  — 

1024  MS 

L  =  256  ( — ]  256  i  Hz  . 

*  W  2 

(A  sample  program  utilizing  (9),  (10),  (11),  (18),  and  the  method  of  appendix  A 
is  given  in  appendix  C  for  75  percent  overlap. ) 

50  PERCENT  OVERLAP 

At  50  percent  overlap  of  the  temporal  weighting,  *  several  possibilities  were 
tried.  They  included 


1 

—  Hz, 
8 


(18) 


cosine  lobe  :  w(t)  =  cos  (»t/lw) 
cosine2  lobe  (Hanning)  ;  w(t)  =  cos2  (jt/Lw) 
Dolph-Chebyshev  (Reference  o)t 


|t|  <  Lw/2  . 


(19) 


A  complete  list  of  cases  is  presented  in  table  1. 

In  figures  4A  through  41,  $  decibel  plots  of  the  magnitude  of  thi  estimated 
spectrum  are  given  for  cosine  temporal  weighting  and  for  (cosine)2  delay 
weighting.  All  plots  are  normalized  with  respect  to  a  maximum  of  0  dB,  which 
occurs  for  f^  =  f  =  256  Hz,  9-  0  Hz.  Figure  4A,  for  example,  demonstrates 
the  behavior  predicted  by  figure  3,  ~amely,  the  presence  of  spurious  sidelobes 
every  S"1  =  (.5xl  s)~l  =  2  Hz.  The  largest  spurious  lobe  in  figure  4A  is 
-23. 5  dB  at  258  Hz.  The  slow  rate  of  decay  of  the  peaks  at  256  +  2n  Hz  is  due  to 
the  discontinuity  of  slope  of  w(t)  at  +  Lw/2  for  this  example.  The  desirable 
feature  of  a  narrow  mainlobe  is  attained,  as  indicated  in  figure  4.  The  succes¬ 
sion  of  plots  in  figure  4  shows  that  the  extent  of  the  picket  fence  varies  greatly 

♦When  these  weightings  are  employed  in  the  FFT,  they  are  delayed  by  1^/2 
seconds,  thereby  being  nonzero  in  the  interval  (0,  L^). 

Ta  quick  and  accurate  method  of  generating  the  Dolph-Chebyshev  weights  by 
means  of  efficient  use  of  an  FFT  is  presented  in  reference  6. 

^Figures  4  through  14  follow  the  text,  beginning  on  page  16. 
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Table  1.  Examples  of  Temporal  and  Delay  Weightings 


Figure 

Temporal 

Weighting 

Delay  Weighting 

Overlap 

(%) 

Number  of  Tones 

4 

cosine 

cosine2 

50 

1 

5 

cosine2 

cosine2 

50 

1 

6 

Dolph-Chebyshev 

cosine2 

50 

1 

I 

cosine2 

cosine2 

75 

1 

cosine2 

cosine2 

75 

1 

4 

cosine 

cosine2 

75 

1 

10 

cosine5 

cosine2 

75 

1 

11 

Dolph-Chebyshev 

cosine2 

75 

1 

12 

cosine2 

flat 

75 

1 

13 

cosine2 

cosine2 

75 

2 

14 

cosine2 

flat 

75 

2 

with  excitation  frequency,  reaching  a  maximum  of  -3. 20  dB  in  figure  4H  for 
=  256  7/16  Hz.  (The  figures  for  %  >  256  1/2  Hz  repeat  the  behavior  shown. ) 
The  worst  sidelobe  of  -23.0  dB  occurs  for  =  256  l/8  Hz,  as  shown  in  figure  4C. 

It  should  be  noted  that  if  sidelobes  were  to  be  measured  with  respect  to  the 
peakonthat  same  plot,  figure  4C  would  yield  a  sidelobe  of  -23.0  +  0. 13  =  -22.9db. 
Thus,  the  convention  adopted  here  must  be  kept  in  mind  in  the  following  discussion. 

Instead  of  applying  the  weighting  directly  in  the  time  domain  by  means  of 
multiplication  on  the  data  x,  the  effect  of  cosine  weighting  can  be  accomplished 
in  the  frequency  domain  by  means  of  convolution  of  the  spectrum  with  the  se¬ 
quence  (i/2)  jl,  -l} ;  however,  the  resultant  must  be  interpreted  as  the  spectral 
value  between  the  two  quantities  convolved  at  each  frequency  step.  *  More  gen¬ 
erally,  (cosine)11  time  weighting  can  be  accomplished  alternatively  by  means  of 
convolution  of  the  (unweighted)  spectrum  with  the  sequence 

©“MS’-® . <-*>1 


♦This  possibility  and  its  interpretation  were  pointed  out  by  Dr.  N.  L. 
Owsley. 
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of  length  n  +  1,  and  then  interpreted  as  the  spectral  value  at  the  center  of  the 
region  convolved,  for  each  frequency  step;  see  appendix  D.  The  convolutional 
sequences  in  (20)  are  given  in  table  2  for  n  =  1  through  5. 

Table  2.  Convolution  .Sequences 


n 

Convolution  Sequence 

1 

i  1/2  {  1,  -l( 

2 

1/4  \  -1,  2,  -1} 

3 

i  1/8  {  -1,  3,  -3,  l| 

4 

1/16 1  1,  -4,  6,  -4,  1| 

5 

i  1/32  j  1,  -5,  10,  -10,  5,  -l| 

Since  the  effect  of  (cosine)2  temporal  weighting  is  very  easy  to  incorporate 
in  the  frequency  domain  by  means  of  convolution,  it  must  be  considered  as  a 
candidate  for  weighting.  The  results  in  figures  5A  and  5B  show  that  although 
the  picket  fence  is  reduced  to  -2.  70  dB,  the  peak  sidelobe  increases  to  -15. 4  dB. 
(For  brevity,  we  are  now  presenting  only  selected  cases  of  worst  excitation 
frequencies. )  The  reason  for  the  increased  eidclebcs  for  this  temporal  weight¬ 
ing  is  that  50  percent  overlap  is  not  yet  great  enough  to  realize  the  deeper  first 
sidelobe  level  of  -31. 5  dB;  that  is,  we  are  still  sampling,  according  to  figure  3, 
on  the  skirts  of  the  mainlobe  for  some  excitation  frequencies.  Generally,  for 
50  percent  overlap,  the  peak  sidelobe  will  occur  approximately  at  the  excitation 
frequency  such  that  the  worst  sidelobe  (or  mainlobe)  of  the  temporal  window 
beyond  f  =  1.  5/Lyy  is  encountered;  this  may  be  seen  by  considering  figure  3 
and  recalling  that  we  plot  only  the  central  portion  of  Y(f, »).  Thus  the  (cosine)11 
weightings  in  tables  1  and  2  for  n  >  2  are  not  acceptable  for  50  percent  overlap, 
since  sampling  of  the  mainlobe  is  encountered. 

The  realization  of  minimum  sidelobe  level  for  a  specified  beamwidth(tothe 
first  null)  is  exactly  the  problem  addressed  by  Dolph,  reference  5.  Accordingly, 
this  weighting  is  of  considerable  importance  in  spectral  estimation.  In  figures  6A 
6C,  the  effects  of  Dolph-Chebyshev  time  weighting  are  presented.  The  worst 
sidelobe  of  -33. 2  dB  occurs  for  f0  =  256  1/8  Hz  (figure  6C).  These  results  are 
noticeably  better  than  in  figures  4  and  5. 

When  triangular  temporal  weighting  was  tried,  it  had  a  peak  sidelobe  of 
-20. 2  dB;  again,  we  are  sampling  the  skirts  of  the  mainlobe.  Thus,  if  50  percent 
overlap  is  all  that  can  be  utilized  for  some  applications,  due  perhaps  to  limited 
computation  time,  the  cosine-lobe  temporal  weighting  is  the  best  of  the  simply 
applied  windows  (that  is,  by  means  of  frequency  domain  convolution),  but  the 
Dolph-Chebyshev  time  weighting  is  10  dB  better  than  the  cosine  lobe  weighting. 
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75  PERCENT  OVERLAP 


At  75  percent  overlap  of  the  temporal  weightings  the  following  examples 
were  utilized: 


cosine2  lobe  (Hanning)  :  w(t)  =  cos2  (#t/Lw) 
cosine2  lobe  :  w(t)  -  cos2  (*t/LyV) 

cosine4  lobe 
cosine®  lobe 
Dolph-  Chebyshev 


:  w(t)  -  cos4  (f  t/lyy) 
:  w(t)  =  cos5  (tt/Lw) 


|t|  <  V2  • 


(21) 


The  results  for  the  Hanning  weighting  are  given  in  figure  7.  The  peak  sidelobe 
is  -41. 8  dB  at  fo  =  256  1/2  Hz  in  figure  7C,  and  the  picket  fence  is  -2. 60  dB 
at  fo  =  256  7/16  Hz  in  figure  7B.  Thus,  a  much  improved  sidelobe  level  is 
realized  relative  to  GO  percent  overlap,  at  the  expense  of  increased  computa¬ 
tion  effort,  that  is,  increased  overlap  and  number  of  FFTs. 


In  an  effort  to  further  improve  performance,  the  (cosine)2  weighting,  which 
has  a  higher  degree  of  continuous  derivatives  at  +  Lw/2,  was  tried.  The  re¬ 
sults  in  figure  8  show  a  maximum  sidelobe  of  -51.  2  dB  and  a  picket  fence  of 
-2. 34  dB. 


Continuation  of  this  effort  to  smoother  weightings  such  as  (cosine)4,  figure  9, 
shows  a  peak  sidelobe  of  -48. 1  dB  and  a  picket  fence  of  -2. 18  dB.  The  peak 
sidelobe  has  increased  over  that  for  (cosine)2  weighting  because,  for  75  percent 
overlap,  the  peak  sidelobe  (or  main  lobe)  of  the  temporal  windowbeyond  approxi¬ 
mately  f  =  3. 5/liyy  is  encountered.  It  so  happens  that  the  worst  case  in  this 
region  is  larger  for  (cosine)4  than  for  (cosine)2  temporal  weighting. 

For  (cosine)5  temporal  weighting,  the  peak  sidelobe  is  reduced  further  to 
-54. 1  dB.  Also,  the  picket  fencs  is  improved  to  -2. 07  dB;  see  figure  10.  This 
weighting  is  easily  accomplished  via  frequency  domain  manipulations;  seetable  2. 
An  alternative  temporal  weighting  of  virtually  equal  quality  to  (cosine)5  is  cubic, 
that  is,  sections  of  cubic  curves  that  have  continuous  derivatives  of  as  high 
order  as  possible.  Ifce  temporal  window  is  proportional  to 


and  has  a  worst  sidelobe  of  -53. 1  dB.  (The  picket  fence  was  not  computed. ) 
However,  a  cubic  temporal  weighting  is  not  easily  accomplished  in  the  frequency 
domain. 


13 


TR  4767 


When  (cosine)6  temporal  weighting  is  considered,  it  is  found  that  75  per¬ 
cent  overlap  forces  us  to  sample  the  temporal  window  on  the  skirts  of  its  main- 
lobe.  This  is  an  unacceptable  weighting  because  the  peak  sidelobe  in  the  vernier 
spectrum  increases  significantly. 

The  possibilities  of  Dolph-Ohebyshev  weighting  are  indicated  in  figure  11. 
The  worst  sidelobe  is  -86  dB  in  figure  11B  for  ?o  =  256  1/4  Hz,  and  the  picket 
fence  is  -2.  29  dB  in  figure  11C  for  fo  =  256  7/16  Hz.  This  is  a  3 2-dB improve¬ 
ment  in  sidelobe  level  compared  with  (cosine)6  weighting.  The  picket  fence  is 
degraded  by  0. 22  dB. 

To  determine  the  effect  of  not  using  delay  weighting,  figure  12  was  computed 
for  (cosine)2  temporal  weighting  and  flat  delay  weighting.  For  certain  excitation 
frequencies,  a  very  narrow  mainlobe  is  realized;  see  figure  12A.  However,  for 
other  excitation  frequencies,  the  lack  of  delay  weighting  creates  broad  "shoulders" 
of  significant  level;  see  figure  12B.  Also,  the  sampling  in  f  and  v,  inherent 
in  the  FFT,  produces  a  picket  fence  of  -5. 0  dB.  Thus,  although  a  low  peak 
sidelobe  is  attained,  lack  of  delay  weighting  is  very  detrimental  to  performance, 
as  will  be  further  demonstrated  below.  Notice  from  figure  7  that  (cosine)2  delay 
weighting  also  yields  a  peak  sidelobe  of  -41. 8  dB,  but  has  no  broad  shoulders 
and  has  a  picket  fence  of  only  -2. 60  dB. 

The  detrimental  effects  of  no  delay  weighting  are  best  illustrated  by  a  com¬ 
parison  of  figures  13  and  14,  which  have  two  tones  separated  by  1/2  Hz,  one 
15  dB  stronger  than  the  other.  It  is  seen  that  these  two  tones  are  resolved,  even 
though  they  are  closer  than  the  resolution  capability  of  the  individual  time  seg¬ 
ments,  that  is,  closer  than  1  Hz.  In  figure  13,  (cosine)2  delay  weighting  is 
employed;  in  figure  14,  none  is  employed.  A  comparison  of  part  A  of  these 
figures  reveals  that,  for  excitation  frequencies  256  and  256  1/2  Hz,  the  flat 
delay  weighting  is  better.  However,  for  excitation  frequencies  256  1/16  and 
256  9/16  Hz,  the  presence  of  the  weaker  tone  is  clearly  evident  in  figure  13B, 
but  hardly  discernible  in  figure  14B  (no  delay  weighting).  The  presence  of  noise 
would  obscure  the  weaker  peak.  Thus,  although  the  peak  sidelobe  may  be  very 
small,  the  presence  of  high-level  broad  shoulders  must  also  be  eliminated  by 
use  of  delay  weighting. 

CONCLUSIONS 

An  approximate  and  quick  FFT  technique  for  vernier  spectral  analysis  is 
possible  by  employing  overlapped  temporal  weighting  and  delay  weighting.  For 
50  percent  overlap  and  Hanning  delay  weighting,  the  best  simply-applied  tem¬ 
poral  weighting  discovered  was  a  single  cosine  lobe,  realizing  a  peak  sidelobe 
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of  -23  dB.  However,  Dolph-  Cheby shev  temporal  weighting  achieves  -33  dB  side- 
lobes.  For  75  percent  overlap  and  Hanning  delay  weighting,  the  best  simply  ap¬ 
plied  temporal  weighting  discovered  was  (cosine)^,  which  realized  apeak  sidelobe 
of  -54  dB.  However,  Dolph- Cheby  shev  temporal  weighting  is  capable  of  -86  dB 
sldelbbes.  Which  overlap  and  weighting  to  employ  depends  on  the  limitations  on 
computation  time  and  storage,  and  on  the  relative  strength  and  location  of  inter¬ 
fering  tones. 

Hie  overlap  is  not  limited  to  the  above  choices.  It  could,  for  example,  be 
67  percent.  Hie  best  weightings  were  not  investigated  in  this  case,  because  it 
was  felt  that  the  above  overlaps  were  easier  to  implement  in  most  cases  of 
practical  interest  However,  Dolph-Chebyshev  weighting  1b  always  a  strong 
candidate  and  is  quickly  and  accurately  generated  (reference  6). 
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Frequency  (Hz) 

Figure  4B.  f0  =  256  1/16  Hz 

4.  Vernier  Spectrum  for  Cosine  Temporal  Weighting,  50  percent  Overlap 
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Frequency  (Hz) 

Figure  4D.  fQ  =  256  3/16  Hz 

Figure  4  (Cont'd).  Vernier  Spectrum  for  Cosine  Temporal  Weighting, 

50  percent  Overlap 
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Figure  4F.  fQ  =  256  5/16  Hz 


Figure  4  (Cont'd).  Vernier  Spectrum  for  Cosine  Temporal  Weighting, 

50  percent  Overlap 
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Figure  4G.  f0  =  256  6/16  Hz 
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Figure  4H.  fo  =  256  7/16  Hz 

Figure  4  (Cont’d).  Vernier  Spectrum  for  Cosine  Temporal  Weighting, 

50  percent  Overlap 
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Frequency  (Hz) 
Fig  re  5A.  fQ  =  256  Hz 
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Figure  5B.  fQ  =  256  7/16  Hz 

Figure  5.  Vernier  Spectrum  for  (Cosine)2  Temporal  Weighting, 

50  percent  Overlap 
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Figure  6A.  fn  =  256  Hz 
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Figure  6B.  f0  =  256  1/8  Hz 

Figure  6.  Vernier  Spectrum  for  Dolph-Chebyshev  Temporal  Weighting 

50  percent  Overlap 
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Figure  6  (Cont'd).  Vernier  Spectrum  for  Dolph-Chebyshev  Temporal 
Weighting,  50  percent  Overlap 
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Figure  7.  Vernier  Spectrum  for  (Cosine)2  Temporal  Weighting,  75  percent  Overlap 
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figure  7C.  ^  =  256  1/2  Hz 

Figure  7  (Cont'd).  Vernier  Spectrum  for  (Cosine)2  Temporal  Weighting, 

75  percent  Overlap 
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Figure  8A.  ^  =  256  Hz 
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Figure  8B.  f0  =  256  1/8  Hz 

Figure  8.  Vernier  Spectrum  for  (Cosine)3  Temporal  Weighting,  75  percent  Overlap 


26 


TR  4767 


Frequency  (Hz) 

Figure  8C.  fQ  «  256  7/16  Hz 

Figure  8  (Cont'd).  Vernier  Spectrum  for  (Cosine)3  Temporal  Weighting, 

75  percent  Overlap 
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246  256  266 


Frequency  (Hz) 

Figure  9C.  ^  =  256  1/2  Hz 

Figure  9  (Coat'd).  Vernier  Spectrum  for  (Cosine)4  Temporal  Weighting, 

75  percent  Overlap 
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Figure  IOC.  ^  =  256  7/16  Hz 

Figure  10  (Cont'd).  Vernier  Spectrum  for  (Cosine)5  Temporal  Weighting, 

75  percent  Overlap 
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Figure  12A.  f0  =  256  1/2  Hz 
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Figure  12B.  f  =  256  7/16  Hz 
Figure  12.  Vernier  plectrum  for  Flat  Delay  Weighting 
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Figure  14B.  f0  =  256  1/16  Hz 


Figure  14.  Vernier  Spectrum  for  Flat  Delay  Weighting,  Two  Tones 
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Appendix  A 


no  +  N- 1 

TWO  METHODS  OF  COMPUTING  ^  exp  (-i2*pn/N)  qn 


n  =  n_ 


Define 


no+N-1 

Op  =  ^  exp  (-i2»pn/N)qn,  0  <p<N-  1  , 

n  “  ^ 

where  >  0.  If  we  let  m  =  n  -  i^,  (A-l)  becomes 

N-  1 

Op  =  exp  (-i2*pn0/N)  ^  exp  (-i2»pm/N)  qm+nQ 

m-0 

.no+N-1 

The  sum  in  (A-2)  is  an  FFT  of  the  sequence  (q^^ 


(A-l) 


(A-2) 


For  an  alternative  method,  consider  the  general  term  qn  in  (A-l).  Then, 


if 


(a)  n  =  0,  N,  2N . qn  gets  weight  exp  (-i0) ; 

(b)  n  »  1,  N+ 1,  2N+ 1, . . . ,  qn  gets  weight  exp  (-i2*p/N)  ; 

(c)  n  =  N- 1,  2N - 1,  3N- 1, . . . ,  qn  gets  weight  exp  (-i2irp(N -  1)/N)  . 
So,  if  we  define  n  =  n  nod  N,  then  case 

(a)  corresponds  to  n  =  0  ; 

(b)  corresponds  to  n  =  1  ; 


(A-3) 


(A-4) 


Therefore,  let 


(c)  corresponds  to  n  =  N  -  1 


vfi=  V  noinino+N" 


(A-5) 


A-l 
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in  which  case 


N  -  1 

v£ 


exp  (-i2t  pn/N)  v~  ,  0<p<N-l; 


(A-6) 


N-l 

this  is  simply  an  FFT  of  |vft  |Q  ,  with  no  phase  factor  necessary.  Equation 
(A-5)  corresponds  to  filling  up  the  array,  from  the  given  quantities  qn, 
starting  from  the  nonzero  position,  iIq  mod  N,  and  cycling  around  to  position  0. 


To  apply  these  results  to  (10),  suppose  weight  w  is  nonzero  for  t  >  0  . 
Then  if  iiq  is  the  smallest  integer  such  that  %  >  k  S/A ,  (10)  canbe  expressed  as 


a 


no  +  N-1 

=  A  ^  ^  exp  (-i2*pn/N)  x(nA)  w(nA  -  kS)  . 
n  =  no 


(A-7) 


This  is  of  the  form  of  (A-l)  if  we  identify 

qn  =  x(nA)  w(nA  -  kS)  . 


(A-8) 


A-2 
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Appendix  B 

DERIVATION  OF  VERNIER  SPECTRUM 
From  the  first  line  of  (7)  in  the  main  text,  there  follows  immediately 

Y(f,  x)  =  A(f, x)  •  D(x)  •  ,  (B-l) 

where  all  convolutions  are  on  v,  for  f  fixed.  Then  using  (6),  we  obtain 

A(f,x)  =y"dr  exp  (-12*  xr)  a(f,r) 

=  / dr  exp  (-12*  vr)  f  dt  exp  (-i2*ft)  x(t)  w(t  -  r)  a$(t)  (B-2) 

J  J  A 

=  W(-x)  J  dt  exp  (-i2*(f  +  x)  t)  x(t)  A{^(t) 

-  W(-»)  -  x)  • 

Substituting  (B-2)  in  (B-l),  we  have 

Y(f,*)=  jw(-x)£x^f  +  x  - 


m 

A 


)} 


►  D(x)  •  5,,(x) 


(B-3) 


B-l/B-2 
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Appendix  C 
SAMPLE  PROGRAM 


The  program  furnished  below  is  illustrative  of  tbe  vernier  technique.  It 
has  been  written  for 


1  1 

A  (DEL)  = - seconds,  N  =  1024,  M  =  32,  S  =— sec 

1024 

fo(F0)  =  256(1/16)  256-i  , 


(C-l) 


but  could  be  ''rsily  changed  to  other  cases.  The  input  data  are  furnishedby  inter¬ 
nal  functions  XREAL  and  XIMAG;  currently,  two  tones  of  relative  strength 
-15  dB  and  separation  1/2  Hz  are  incorporated.  Loop  1  in  the  main  program 
accomplishes  Hanning  temporal  weighting,  while  loop  2  accomplishes  Hanning 
delay  weighting.  The  subroutines  MKLFFT  and  QTRCOS  are  detailed  in  refer¬ 
ence  7.  The  method  in  this  program  employed  the  cycling  technique  described 
in  appendix  A. 


parameter  n=i oi<*  * m=32 ,  n4=n/4*i  * m*sm/4*i 
integer  ps 

DIMENSION  ZR(N).ZI(N),«(N)»D(M)»AR121,M)*AI(21*M)»ADR(M)*ADKM)* 
SDQ(M)  »CN(N4) »CM(M4)  ,Z(200)  rX(it>6)  »Y(168) 

S8=l0»**(-.75> 

7PIMS2 , *3,141592654/ ( M+l ) 

TpIN=t,*3,l41592654/N 

M1=M-1 

N1SN-1 

DELsl./N 

ISSN/* 

SslS*PEL 

IMSINTIL0G( FLOAT (M) >*1.4427*. 5) 

CALL  8TRC0S(CN»m) 

CALL  8TRC0S(CM»v|) 

call  modess(2»o j 

call  suaJEG(z»a.»-ioo,,i65,,o,) 

CALL  0djCTS(2#ea0t, 1400, »2425,, 2350.) 

DO  1  MSsOtNl 

1  W(MSU)s1,-C0S(tPIN*MS) 

DO  2  KS=0»Ml 
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2  D(KS+1)=1.-C0S(tPIM*(KS+1>) 

Do  13  IC=lrl6B 

13  X(lC)fZC 

Do  12  IF=0.8 

F0S256.+IF*,0625 

00  3  KSsOrMl 

N0=IS*KS 

Nu=N0+N1 

00  4  NSsNOfNU 

NTSMOO(NS»N) 

ZR ( NT+1 ) =XRE AL ( NS* jEL ) *  W ( NS-NO+1 ) 

Zl  (NT+i )  =X1MA6  (  nS*u£L  )  *N  ( NS-NO+l ) 

4  CONTINUE 

CAUL  HKLFFT (ZR*Zlf tNf 10»"1) 

00  5  PS=246,266 

AR (PS-245* KS+1 >s  DtL*ZR(PS+l ) 

AI(pS-245.KS+l)=  OtL*ZI (PS+1 ) 

5  CONTINUE 

3  CONTINUE 

00  6  PS=24&,266 
00  7  KSsOtMl 

ADR (  KS+1 ) =&• AR ( pS-245 , KS+1 ) *0  (  KS+1 ) 

AOI (  KS+1 » =S*AI ( pS-245 , KS+1 ) *0 ( KS+1 ) 

7  continue 

CALL  MKLFFT(ADR.A0I,CM»IM»-1) 

IFUF.EQ.O.ANO.pS.EQ.256)  PEAKs1Q.*L0310(AQR(1)**2+ADK1)**2> 
00  8  KSsO.Ml 

AR ( PS-245*  KS+1 ) = ADR (KS+1 1 
A I ( PS-245  *  KS+1 ) sADI ( KS+1 > 

8  CONTINUE 

6  CONTINUE 

00  9  PSs246# 266 
ICS (PS-246) «8 
00  10  MS=0*M1 

AsAR (PS-245 * HS+i }**2+AI (PS-245*MS+1) **2 
AsMAX (A*l*E-36) 

OB (MS+1 ) slO , *L0si0 ( A) -PEAK 

10  continue 

PRINT  14.  PS 

14  format (///i io/) 

PRINT  llf  08 

11  format  (/Bins. 6) 

Y(IC+ii=0B(M-3) 

Y(1C+2)sOB(M-2) 

Y(1C+3)sOB(M-1) 

Y(IC+4)=0B(M) 

Y(IC+5)sOB(l) 

Y(IC+6)=0B(2) 

Y ( IC+7) =DB(3) 

Y(IC+8)s08(4) 

9  CONTINUE 
SlOEs-200* 

00  17  IC-lt 168 

IF(IC,GE*69, ANO. IC.LE, 109)  80  TO  17 
S|UE=MAX(SIDE»t(IC)) 
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17  CONTINUE 
PRINT  18*  SIOE 

18  FORMAT t/'  PEAK  SIDEL086  IS*,E14.S) 

call  setsm6(z»3o»i.) 

DO  15  IP=10#90» j 0 

call  kiNES6<z»o.  b. ♦ -float cxp> ) 

15  Call  kXNES6<Z* 1,165.* -float cxp* ) 

DO  16  IPsl3,157.8 

call  lxnes6(z»o,floatup)»o«) 

16  call  lxnes6<z»i. float (xp>*-ioo.> 

Call  S£TSM6(Z»3n*2.) 

call  kiN£S«(z*o.5.,-ioo.> 
call  lxness(Z»i,5.»o.) 

Call  kXNES&(Z»1.16b.*0.) 
call  kiNESS(z*i.i6b.»-ioo.) 

Call  kXNESS(Z*1.5.»-l00.) 

CALL  kXNES6(Z*l68»X*Y) 

CALL  PA6E6(Z*0*  3*1) 

12  CONTINUE 

call  Exits(z) 

PRINT  11*  PEAK 

function  xreal(t) 

XREAL»C0S(2.#3.i41S>92654*F0*T) 

S  ♦CLS(2.#3.l41592654*(F0*.5)*T)*S9 

return 

function  ximasit) 

XlMAS»SIN(2.*3.l41t,92C54*F0*T) 

S  ♦SIN(2.#3.i41b92654*(F0t.t»)*T)*SQ 

return 

EnO 


C-3/C-4 
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Appendix  D 

EFFECT  OF  (COSINE)11  TIME  WEIGHTING 

Let  us  define  the  spectrum  of  unweighted  data  x  as 

N-  1 

Z(f)  E  fdt  exp  (—12 » ft)  x(t)  A5A(t)  =  A  £  exp  (-i2*kAf)  x(kA),  (D-l) 


(D-3) 


and  that  corresponding  to  weighting  w  as 

Z^f)  =  fdt  exp  (-12*  ft)  x(t)  w(t)  A6A(t)  =  Z(f)  •  W(f) 

N-l  (°-2) 

=  4  V  exp  (-i2*kAf)  x(kA)  w(kA)  . 
k  =  0 

Now  for  (cosine)11  time  weighting,  we  have* 

w(kA)  =  sin11  (k»/N),  0  <  k  <  N  -  1  •  (D"3 

Substituting  (0-3)  in  (D-2),  there  follows 

n  *1  ^ 

zw(f)  =4  exp  (-12.k4f)  x(k4)  Q  [exp  (Ik./N)  -  exp  (-Ik,/N)J 

-t— «fe)‘ 


■0  fr-o  ■{-=) 


(D-4) 


the  last  step  via  use  of  (D-l).  This  result  yields  (20). 

- *See  the  firstfootnote  to  equation  (19)  of  the  main  text  for  the  explanation 

of  sinn  in  (D-3). 
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For  the  special  case  of  n  =  1,  (D-4)  becomes 


and  in  particular, 


Sw(f>  =  i  1/2  |z(f-  —  )  -  Z[f+  ~  ]|, 

L  \  2NA/  \  2NA/J 

— +  — )=i  1/2  Tz  (~y  z  (— )] . 

NA  2NA/  [  V  NA  /J 


(D-5) 


(I>6) 


The  right-hand  side  of  (D- 6)  involves  two  adjacent  spectral  values  as  afforded 
by  the  standard  N-point  FFT  in  (D-l).  The  left-hand  side  of  (D-6)  is  the  spec¬ 
tral  value  of  2^  at  the  frequency  halfway  between  the  above  two  spectral  loca¬ 
tions.  Thus,  under  this  interpretation  of  the  right-hand  side  of  (D-6),  the  desir¬ 
able  sidelobe  control  predicted  by  the  convolution  in  (D-2)  can  be  attained.  A 
similar  interpretation  is  possible  for  (D-4)  for  general  n. 
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